Effect of local charge fluctuations on spin physics in the Neel state of La 2 Cu04 
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We explore the effect of local charge fluctuations on the spin response of a Mott insulator by 
deriving an effective spin model, and studying it using Schwinger boson mean field theory. Applying 
this to La2CuC>4, we show that an accurate fit to the magnon dispersion relation, measured by Coldea 
et al. [Phys. Rev. Lett. 86, 5377 (2001)] is obtained with Hubbard model parameters U « 2.34eV, 
and t ~ 360meV\ These parameters lead to estimates of the staggered magnetization (m s ~ 0.25), 
spin wave velocity (c ta 800mel/-A), and spin stiffness (p s ta 2AmeV). In particular the staggered 
moment as well as the effective local moment are renormalized to smaller values compared to the 
Heisenberg model due to local charge fluctuations in the Hubbard model. The dynamical structure 
factor shows considerable weight in the continuum along the zone boundary as well as secondary 
peaks that may be observed in high resolution neutron scattering experiments. 

PACS numbers: 



I. INTRODUCTION 

Multispin interactions are important in a variety of 
magnetic systems^ In a magnetically ordered insulator, 
such interactions can reveal themselves through their ef- 
fect on the spin dynamics measured in neutron scattering 
experiments. One system which appears to provide an 
example of such physics is La2Cu04, which is a Mott in- 
sulating antiferromagnet. Indeed, high resolution data on 
the magnon dispersion in this system indicate a dip in the 
magnon energy near Q = (tt/2, 7t/2) while traversing the 
magnetic Brillouin zone boundary along (77, 0)— >(0, tt). 
Such a dip is not expected for the magnon dispersion 
in a nearest-neighbor Heisenberg antiferromagnet, but it 
has been suggested^ that this could arise from four-spin 
interactions generated by local charge fluctuations in the 
insulating phase. Motivated by this, we revisit the zero 
temperature Neel ordered state of La2Cu04, and exam- 
ine how the spin stiffness, the ordered moment, the ef- 
fective local moment at a site, and the spin dynamics 
are influenced by local charge fluctuations in this Mott 
insulator. 

La2CuC>4 is a layered antiferromagnet. While the ex- 
change coupling between the two-dimensional layers is 
crucial for the existence of a finite-temperature Neel tran- 
sition in this system, this interplane coupling is neverthe- 
less tiny in magnitude, ~ 10~ 5 of the in-plane exchange 
coupling, and much smaller than the resolution of the 
neutron scattering experiments. It is therefore expected 
that the low temperature spin dynamics in the Neel or- 
dered state in this material is adequately captured by 
the Hubbard model at zero temperature on a two dimen- 
sional square lattice, (with standard notations) 

H = - t J2 ( c I c j° + h - c ) + U Y1 n ^ w 

in its insulating phase at strong coupling U ^> t. In this 
regime, it is well known 3 ^ that the low energy physics can 



be encapsulated in an effective spin Hamiltonian which 
incorporates two-spin and four-spin interactions (given 
below in Eq. [2J . 

A detailed analysis of a spin model with such four-spin 
interactions was carried out by Katanin and Kampfe us- 
ing modified spin wave theory. This allowed them to con- 
clude that a fit to the magnon dispersion indeed requires 
sizeable four-spin interactions, as originally suggested by 
Coldea et a& based on linear spin wave theory. A differ- 
ent approach to this problem was adopted by Singh and 
GoswamiS and Peres and Araujoi who used the random- 
phase approximation (RPA) for the Hubbard model with 
a spin density wave ground state. This is presumably a 
better approximation at weak coupling but nevertheless 
also reproduces the dip in the magnon dispersion near 
(71-/2, 7r/2) at intermediate values of U/t. 

In this paper, we use Schwinger boson mean field the- 
ory to study the magnon dispersion and scattering con- 
tinuum for the effective spin model in Eq. © . There are 
several reasons to try this route. First, in contrast to the 
RPA studies, this approach of mapping to an effective 
spin model is explicitly geared towards addressing the 
strong coupling limit U/t 3> 1. It is therefore interesting 
to compare the two approximations at the intermediate 
coupling values of experimental interest, especially taking 
into account modifications to the spin operator which re- 
sult when going from the Hubbard electronic description 
to an effective spin model. Indeed, it is satisfying that 
some of the results obtained here are consistent with the 
spin-density wave approach for intermediate U/t. Sec- 
ond, Schwinger boson mean field theory is known to work 
well for the square lattice Heisenberg model and it is 
worth asking if it continues to be reliable once sizeable 
four-spin interactions are included, and how such interac- 
tions affect the ground state properties and the magnon 
dispersion. Finally, we consider the scattering continuum 
beyond the magnon peak in La2Cu04 which has not been 
addressed sufficiently in these earlier papers. (While this 
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paper was being refereed, a closely related workS which 
discusses the sum-rule and other issues with neutron scat- 
tering experiments in La2CuC>4 was submitted for pub- 
lication. We have therefore added a brief discussion of 
sum-rules within our approach, and refer the reader to 
the above paper for a detailed analysis of sum rules, the 
current status of neutron scattering experiments and al- 
ternative theoretical approaches to spin dynamics in this 
system.) 

A summary of our main new results is given below. 

(i) Comparing the mean field theory results for the 
ground state energy and static spin structure factor with 
exact diagonalization studies on small clusters, we show 
that the mean field theory provides a good description of 
the ground state of the spin model J2J , in the experimen- 
tally relevant regime. 

(ii) We show that the experimental magnon disper- 
sion is well reproduced by our approach for Hubbard 
model parameters U ~ 2.34eV, and t as 360meV. These 
parameter values are in reasonable agreement with the 
experimental estimates* of U, and electronic structure 
calculations 10 for t. Our estimated ratio U/t « 6.5 is in 
rough agreement with the value inferred from RPA calcu- 
lations for the Hubbard model^ as well as the quantum 
Monte Carlo (QMC) calculations of Sengupta et aLjii 
though we do not resort to the single mode approxima- 
tion used by the latter authors. It is however slightly 
smaller than the ratio U/t w 8.8 obtained from a recent 
series expansion study of the Hubbard modeli^. 

(iii) Our best fit parameters for t, U lead to values for 
the staggered magnetization m s « 0.25, spin wave veloc- 
ity c w 800 meV-A, and spin stiffness p s w 24meV. The 
spin wave velocity and spin stiffness are in good agree- 
ment with experimental findings^ The staggered mo- 
ment is substantially reduced from that of the Heisen- 
berg model due to local charge fluctuations, as was ear- 
lier pointed out by Delannoy et a/, 14 and it may be worth 
revisiting this in experiments. 

(iv) Turning to the spin excitation spectrum beyond 
the magnon, we find that there is a broad continuum of 
excitations at energy scales relevant to neutron scatter- 
ing experiments. In addition, the spin-flip operator in the 
Hubbard model induces transitions into the upper Hub- 
bard band which would show up at energies of 0(U). 
We however restrict attention in this paper to energies 
well below the Mott gap, where a spin model description 
is appropriate. In this regime, we find that, along the 
zone boundary, the spectral weight in the continuum is 
about 40% of the total spectral weight. This is similar to 
the large continuum spectral weight inferred from QMC 
calculations on the Heisenberg modelji^ but we cannot 
make any quantitative comparisons at this stage. Turn- 
ing to the spectral lincshapc at energy scales relevant to 
neutron scattering experiments, we find that the dynam- 
ical structure factor exhibits secondary peaks within the 
mean field theory. 

(v) The local charge fluctuations in the Hubbard model 
also lead to a smaller effective local moment. To ob- 



tain this effective moment, we consider integrating the 
dynamical structure factor over energies well below the 
Mott gap, as relevant for neutron scattering experiments, 
and then integrating this "low-energy" equal-time struc- 
ture factor <Si ow (q) over q. This leads to the sum-rule 
E q ^iow(q) = S cS (S cS + 1), with S cS < S = 1/2. Our 
calculation gives S c g ~ 0.39. 

It is possible that all these features could be explored 
in neutron scattering experiments with higher resolution 
and intensity in the near future. 



II. EFFECTIVE SPIN MODEL 

The effective spin model describing the strong coupling 
regime of the Hubbard model in takes the form 

#spi„ = \Y. J ^ - r *)& • + Yl Jn [& ■ ■ S 4 ) 



+ (Si ■ S 4 )(& • S 3 ) - (Si • S 3 )(S 2 ■ Si) 



(2) 



Here, Si are spin-half operators, and Jn refers to an el- 
ementary plaquette on the square lattice with 1 — 4 la- 
beling sites on its corners. The explicit expressions for 
the exchange couplings in terms of the Hubbard model 
parameters are 4 ", 
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With hopping matrix elements between further neighbor 
sites in the Hubbard model, we expect additional cou- 
plings in the effective spin Hamiltonian. However, elec- 
tronic structure calculations^ suggest that these matrix 
elements are small in magnitude for La2Cu04, so we will 
not consider them here. 



III. SCHWINGER BOSON MEAN FIELD 
THEORY 

To study the ground state and excitations of -f/spin, we 
follow Ref. 16 and represent spins using two species of 
bosons bi 2 as 



Si = b\ a a a ab 



with the constraint 



K a b ia = 2S 



(7) 



(8) 
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on the boson number at each site. The a are Pauli matri- 
ces, the spin 5=1/2 for our system and, unless specified, 
the summation over repeated Greek indices is implied 
here and below. The Hamiltonian can be reexpressed in 
terms of the bosons by using the identity 



Si • Sj 



B h B a 



-a\ a a. 



ij ij i 



(9) 



where :: is the standard normal ordering and the bond 
operators are 



^3 = ~ b L b ja 



(10) 
(11) 



In the exact Schwinger boson representation, the Hamil- 
tonian obtained using the above identity consists of 4- 
boson and 8-boson operators, and one has in addition 
to deal with the boson number constraint JSJ. In order 
to make progress we can resort to a mean field approach. 
Applied to the Hamiltonian H sp i n in Eq. J2J , this consists 
of three approximations, (i) First, we reduce the 4-spin 
interactions to effective 2-spin terms by setting 



(Si • S 2 )(S 3 ■ Si) 



<(Si •S 2 ))(S 3 -S 4 ) (12) 
(Si • S 2 )((S 3 • S 4 )) (13) 
((Si-S 2 ))((S 3 -S4)). (14) 



Defining 
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= J° ff (r) J B(r)cos(k-r) 
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(18) 




= ^ J cff (r)A(r)e^ lk - r 
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B(k) 
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and choosing 



fik 


= y/\Pk - M 2 - l"k| 2 


(22) 


»k 


= cosh(6l k )e-^ k 


(23) 


i'k 


= sinh(6» k )e- 47k 


(24) 


cosh(26» k ) 


= |/3 k -A|/fi k 


(25) 


sinh(26> k ) 


= |«k|/^k- 


(26) 



we arrive at the diagonal form 

= ^/^/ k ^ + A(S + l/2)7V sito 
- l^/3(k) + (A*(k)a k -B*(k)(3 k ) (27) 



This renormalizes the two-spin exchange couplings J(ji — 
Yj) — > J cff (n — rj) in -ff S pin. It is straightforward to show 
that Jf s — J\ + 2J a {Si-Sj) where ij are nearest neighbor 
spins, Jf = J2 — Jn (Si ■ S m ) where im are next near- 

TeF 
y .3 



est neighbors, and J§ = J%- At this stage, we obtain a 



Heisenberg spin Hamiltonian -ff ™[ n with only two-spin in- 
teractions, (ii) Next, we do a mean field decoupling of the 
4-boson terms arising from the effective two-spin interac- 
tions, setting : B%B^ (B\ ] )B l] +B\ ] (B. ll )-(B\ 1 )(B l] ) 
and similarly for A\-A^. This reduces the spin Hamilto- 
nian if™[ n to a quadratic boson Hamiltonian, which up 
to overall constant terms is given by 



off 

boson 



id 

- (4>4-4(4') 



(15) 



(iii) Finally, we take the constraint JSJ) into account 
using a site-independent Lagrange multiplier A, thus 
working with the Hamiltonian H{^ son — H^ son — 



X T,i (b\ a b ia -2S 
diagonalized by a Bogoliubov rotation 



The Hamiltonian H^ Bon 



6 k,i = M k/ kT - «k/l ki 

-k,2 — _U k/kT "+" U k/-kl' 



is readily 

(16) 
(17) 



The /-particles appearing here are bosonic S = 1/2 
spinous. The above parameters u k , w k depend on the 



mean field values of (fly) = B(r), (Ai 



A(r), and 



(Si ■ Sj) appearing in the quadratic Hamiltonian. These 
mean field parameters are evaluated in the ground state 
of 1|27|) , which results in a self-consistent theory. For com- 
pleteness, the self consistency conditions and the con- 
straint equation are given by 
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B(k) 
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E 



(28) 
(29) 
(30) 



It is well known — that long-range magnetic order appears 
in this formulation if the energy r2 kjM vanishes at some 
wavevector(s) k = {Q;} in the thermodynamic limit, 
leading to condensation of the spinons at these momenta. 
In the Neel state on the square lattice, the spinons con- 
dense at wavevectors Q = ±(ir/2, tt/2) in the thermody- 
namic limit 17 leading to magnetic order at the wavevector 
connecting these points, namely (tt, tt). We report below 
the ground state properties of the model with four-spin 
interactions, and present a comparison with exact diag- 
onalization results on small clusters. 
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FIG. 1: Comparison of Schwinger boson mean field theory and 
exact diagonalization calculations of the ground state energy 
per site for the spin model in Eq. J3J for the Heisenberg limit 
(U/t-* oo) and for U/t = 6.5. 



IV. GROUND STATE PROPERTIES 

A. Comparison with exact diagonalization results 
on small clusters 



a; 
in 




FIG. 2: Comparison of Schwinger boson mean field theory 
and exact diagonalization calculations of the static (equal- 
time) spin structure factor S(q) on a 4 x 4 lattice. These 
results are for the bare spin-spin correlation function, where 
the transformation of the spin operator in going from the 
Hubbard model to the spin model have not been included — 
including them modifies the mean-field and the exact diag- 
onalization results in the same manner, as discussed in the 
text, and does not affect our comparison. 



Schwinger boson mean field theory provides very accu- 
rate results for the ground state properties of the nearest- 
neighbor Heisenberg antiferromagneti£. In order to as- 
sess the accuracy of this approximation scheme in pres- 
ence of the four-spin interaction term it is useful to com- 
pare the mean-field estimates of ground state proper- 
ties with exact diagonalization numerics. To this end, 
we have performed Lanczos exact diagonalizations of the 
spin Hamiltonian © on the 4x4, y/20 x y/20, V32 x V32 
and 6x6 clusters. In particular we have focused on 
exchange interactions corresponding to U/t — ► oo (the 
Heisenberg limit), and U/t — 6.5, for which the mean 
field theory fits the experimental magnon dispersion of 
La2CuC>4 as discussed in the next section. The latter 
choice of U/t corresponds to J2 — >h — 0.0276Ji and 
Ju = 0.55 Ji in 0. 

We have calculated the ground state energy and the 
static spin structure factor S(q) 

S(q)=^<5(r i )-S(r ))e- it >-^- r °) (31) 
3 

The results presented in Figs. H an d El clearly show that 
the Schwinger boson approach is still rather accurate 
even for large four-spin interactions which are generated 
for a moderate interaction strength U/t = 6.5. We next 
turn to the staggered magnetization and spin stiffness in 
the thermodynamic limit for this parameter value, tak- 
ing into account the modification of the spin operator by 
local charge fluctuations. 



B. Staggered magnetization and spin stiffness in 
the thermodynamic limit 

As mentioned earlier, the local charge fluctuations in 
the Hubbard model affect the spin correlations. We can 
take this into account by noticing that the spin operator 
of the Hubbard model is modified by the same unitary 
transformation which enables us to rotate from the Hub- 
bard model Hamiltonian to the Heisenberg model with 
four-spin interactions of Eq. which has no charge fluc- 
tuations. This transformation was discussed earlier in 
the context of various physical correlations in the doped 
cuprate superconductors^, and its importance for spin 
physics in the undoped Mott insulator was pointed out 
by Delannoy et ali^ 

For instance, making the unitary transformation on 
the spin operator S + (r), we find 

S+(r) = 5 + (r)(l-4^) + ^ £ S+(r + a) 

a — ±x,±y 

~ Jj H ( C r +Q , T C ra - C r, T C r+aa) ( 32 ) 
a— ±x,±y 

While the first two terms lead to spin flips within the 
lower Hubbard band, the last term excites an electron 
into the upper Hubbard band. This final term plays a 
role at high energies, of 0(U), and short distances^. We 
will therefore not consider it below, but restrict atten- 
tion to energies well below the Mott gap, where the first 
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two terms alone play a role and our spin-only description 
is valid. With this caveat, we find at 0(t/U) 2 we find 
that the transformed "low energy" equal-time spin-spin 
correlation function is given by 



<Siow(r - r') 



(S(r).S(r'))(l 

f_ 

U 2 



a) ■ 5(r')) 



E 0(r- 

a— ±x,±y 

+(5(r) • S(r' + a)) (33) 



Within Schwinger boson mean field theory, we obtain 
the antiferromagnetic order parameter m s from finite-size 
scaling of the static structure factor, defined via 



5iow(q) = y^<Siow(rj - r )e 



as 



lim 

Af sito ^oo 



<5low(7T,7r) 



site 



(34) 



(35) 



We find that the staggered magnetization for U/t = 6.5 
extrapolates to m s w 0.25, smaller than that for the 
nearest- neighbor Heisenberg model (m" n ~ 0.303). This 
decrease of the staggered magnetization is due to lo- 
cal charge fluctuations in the insulator, and arises from 
0(t/U) 2 terms in the modified spin operator in Eq. i|32|) 
as pointed out in Ref. IbH 

For a pure spin model of the copper sites in La2Cu04, 
we expect the sum-rule on the equal-time structure fac- 
tor E q S(q) = s ( s + !) with s = V 2 - However, with 
the above modification to the spin correlations, this sum- 
rule is corrected by terms of order (t/U) 2 . Specifically, 
imagine integrating the dynamical structure factor over 
energies relevant to neutron scattering experiments, i.e. 
much smaller than the Mott gap, and then integrat- 
ing this "low-energy" equal-time structure factor iSi ow (q) 
over all q. In this case, we find the modified relation 
E q Siow(q) = S eS {S eS + 1) with S cff ~ 0.39. It would 
be interesting to look for this apparent sum-rule deficit 
in neutron scattering studies. 

Following Stringari— , we extract the spin stiffness from 
the small momentum behavior of the structure factor. 
Specifically, approaching the magnetically ordered state 
in a rotationally invariant formulation, the static struc- 
ture factor for small q is related to the spin stiffness 
through 



Ps = 



c5(q 0) 



q 



(36) 



where c is the spin wave velocity of the linearly dispersing 
magnon at momenta near (0,0) and (7r,7r). The spin- 
wave velocity c can be obtained from the magnon dis- 
persion as discussed in connection with the excitation 
spectrum in the next section. This, together with the 
calculated structure factor in the limit, q — > 0, yields the 



spin stiffness through Eq. I|36|l . As discussed in detail in 
the next section, a fit to the overall magnon dispersion 
in La2CuC>4 yields a spin wave velocity c ~ 800meV-A. 
Assuming a lattice spacing 3.85A for the CuC>2 plane, 
this yields p s ~ 2AmeV , which is consistent with experi- 
mental findings^ Note that the spin stiffness can also be 
obtained by studying the change in ground state energy 
for a slowly varying spin twist. 



V. EXCITATION SPECTRUM 

To probe the spectrum of S = 1 excitations at a 
given momentum q, which is of interest for neutron 
scattering experiments, one can act with the operator 

5' + (q) = Ek & iL A-q,2 onto the ground state22i22. Mak- 
ing a Bogoliubov rotation to spinon variables, we can 
rewrite this as 



k 

X ("k-q4_ q , T - Uk-qjl k+q , J 



(37) 



Including corrections to the spin operator via the uni- 
tary transformation, discussed earlier for the static struc- 
ture factor, leads to two modifications in the dynami- 
cal response. The first is a simple modification of the 
spin operator by a multiplicative "form factor" S + (q) — > 
G(q)S' + (q), with 

2t 2 

G(q) = 1 - jjj (2 - cos q x - cos q y ) (38) 

This will not affect the magnon energies but will cor- 
rect the dynamical response by an overall multiplicative 
prefactor G 2 (q). 

The second effect of this unitary tranformation as dis- 
cussed is to generate a new term which causes a transition 
from the lower to the upper Hubbard band. This term 
plays a role only if we examine spin dynamics at very 
high energies, and is of no relevance at energies probed 
in the neutron scattering experiments to which we re- 
strict attention here. We see that at a given momentum 
q, if one of the two spinons combining to give the spin 
operator in Eq. (|37() is condensed (which would happen 
for spinon momenta ±(7r/2, ir/2) in the Neel state), that 
part behaves as a single particle excitation with a well de- 
fined dispersion - this is the magnon. For general k, both 
spinons are uncondensed and this remainder of the sum 
contributes to the scattering continuum. Quite gener- 
ally, both parts play a role when we evaluate the spectral 
function for the £ + (q) operator. 



A. Magnon dispersion 

From the above discussion, it is clear that the magnon 
energy at momentum q within the mean field theory 
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FIG. 3: Fit of the magnon dispersion obtained from Schwinger 
boson mean field theory to the experimental data at T = 
WK in Ref. along the standard contour in the Brillouin 
zone. The Hubbard model parameters used for this fit are t = 
360meV and U = 2.34eV , with the effective spin Hamiltonian 
couplings given by Eq. JSJ. 



of the Neel state is just given by r2q_ q /2 where Q = 
±(tt/2, 7r/2). For our choice of mean field decoupling, 
and for spin S — > oo, this is exactly half the value given 
by spin-wave theorjiSi. It is known2£ that Schwinger bo- 
son mean field theory for the Heisenberg model has this 
shortcoming, which can be fixed by working within a 
large- N generalization of the Schwinger boson theory and 
including 1/N corrections to the mean field (N = oo) re- 
sult. Since we reduce our original Hamiltonian in Eq. (2) 
to an effective Heisenberg-type model before we use the 
Schwinger boson representation, we expect the same fluc- 
tuation corrections to appear in our case as well, however 
an explicit calculation is beyond the scope of this paper. 
Here and below we will work with f^k as the magnon en- 
ergy. Note that the 0{t/U) 2 corrections to the spin op- 
erator do not affect the magnon dispersion. Fig. [3] shows 
the magnon energy fl^ obtained within the mean field 
theory for U/t = 6.5 and t = 360meV compared with 
the experimental data from Coldea et al. This choice of 
parameters gives a good fit to the dispersion over the 
entire contour along the Brillouin zone, including both 
the linearly dispersing spin-wave regime as well as the 
high energy magnon dispersion along the zone boundary. 
While the fit along the zone boundary is easy to see in the 
figure, its accuracy at low energies may be seen from the 
spin wave velocity in the mean field theory c ~ 800mey- 
A being in good agreement with experiment. We also 
point out that at much larger values of U/t 7), we are 
unable to reproduce properly the dip in the magnon dis- 
persion at (7r/2, 7t/2) while traversing the zone boundary. 
While further neighbor hoppings in the Hubbard model 
(e.g. third-nearest neighbor hopping £3) could produce 



a similar dip near (7r/2,7r/2) by generating an antifer- 
romagnetic exchange of order t^/U, electronic structure 
calculations for La2Cu04 indicate that such hopping ma- 
trix elements arc too small to be relevant. 



B. Scattering continuum 

Having fixed the parameter values U, t from fitting the 
magnon dispersion to the data, we next turn to the scat- 
tering continuum. Specifically, we are interested in ask- 
ing how much weight is present in the continuum relative 
to the magnon, and if there are any features in the con- 
tinuum which may be experimentally observable. Recall 
that in the mean field theory this part of the spectrum 
is just the two-spinon excitation continuum. In order to 
obtain the weight in the continuum versus the magnon, 
we consider the energy integrated response. This is just 
the structure factor 



S(q) = £ 



|/?k-q - A| 



k 

«k a k-q 



n 



1 



k-q 



1 



k"k-q 



(39) 



To obtain the magnon weight part, we need to keep only 
those contributions where fik or f2k-q would vanish in 
the thermodynamic limit (or are minimum on a finite 
but large lattice). The continuum contribution is the re- 
maining part of the sum. Fig. 0] shows the ratio of these 
weights along the same contour in the Brillouin zone over 
which the magnon dispersion is displayed. This ratio is 
clearly insensitive to the "form factor" G(q), however the 
important thing to note from the plot is that the contin- 
uum weight is most significant along the zone boundary 
and accounts for nearly 40% of the total spectral weight. 
This is the most promising region to study the contin- 
uum in neutron scattering experiments. The continuum 
weight is also considerable along the line from (ir, 0) to 
(■7T, 7r/2). This is qualitatively similar to results inferred 
for the Heisenberg model from QMC calculations^. In 
order to display possible interesting features in the con- 
tinuum scattering we have plotted the spectrum of ex- 
citations at a few q points along the zone boundary in 
Fig. El as well a grayscale plot of the spectral function 
in Fig. O along (0,0) — * (7r,7r) and along (n, 0) — ► (0,7r). 
We see that in addition to the magnon, there are sin- 
gular secondary features in the spectrum. These sec- 
ondary peaks can be shown to arise from special points 
{K} in the Brillouin zone in the vicinity of which the 
sum fik + fik-q varies slowly (dispersing as ~ |k — K| 4 ) 
giving rise to a log singular density of states for the two- 
spinon continuum. These secondary peaks and contin- 
uum scattering also arise within the mean field theory 
for the nearest neighbor Heisenberg model, with minor 
quantitative changes due to differences in the spinon dis- 
persion. The effect of fluctuations beyond the mean field 
result on the magnon dispersion and the scattering con- 
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FIG. 4: Ratio of continuum and magnon spectral weights 
along a contour in the Brillouin zone. The magnon exhausts 
the spectral weight for q — » (0, 0) and q — » (n, tt) consistent 
with general arguments—. Along the zone boundary, the con- 
tinuum accounts for nearly 40% of the total spectral weight. 



q= (71/4,371/4) 
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FIG. 5: Grayscale plot of the dynamical structure factor 
S(q, E) of the spin model with optimal parameters as in 
Fig. 1, with high intensity in black. The gray areas indicate 
regions of continuum scattering. The intensity has been scaled 
to clearly show both the magnon and secondary peak features. 
Left panel: Magnon peak and the single secondary peak along 
(0, 0) — > (n, 7r). The magnon and secondary peak intensities 
vanish as q — > (0, 0). Right panel: Magnon peak and two sec- 
ondary peaks along the zone boundary (7r,0) — > (0, tt). The 
two secondary peaks merge at (0,7r) and (^,0). The dip in 
the magnon dispersion at {tt/2,-k/2) is clearly visible. 
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FIG. 6: The dynamical structure factor along the zone bound- 
ary at (7r/2,7r/2) (top), (tt/4, 3tt/4) (middle) and (tt, 0) (bot- 
tom), showing the sharp magnon (as a bold vertical line not 
drawn to scale) and the broad continuum including the sharp 
secondary features most clearly visible at (tt/4, 37r/4). 



tinuum is beyond the scope of this paper, and is currently 
being investigated. 



VI. CONCLUDING REMARKS 

We have presented a Schwinger boson mean-field de- 
scription of the of the effect of four-spin interactions on 
spin dynamics in Mott insulators. Such multispin inter- 
actions and their effects are potentially important in any 
correlated insulator with local charge fluctuations. For 
La2CuC>4, the experimental magnon dispersion is well re- 



produced by our approach with Hubbard model parame- 
ters U « 2.34eV, and t ps 360meV. These values are con- 
sistent with the experimental estimates^, and electronic 
structure calculations^. This leads to a sizeable value of 
the four-spin exchange term J a ~ 0.55Ji in the effective 
spin Hamiltonian. Good agreement with experiments is 
also obtained for spin wave velocity and spin stiffness, 
thus leading to a consistent description of the magnetic 
behavior of this Mott insulator. The staggered moment 
is substantially reduced from that for the Heisenberg 
model, it would be worth revisiting experiments to test 
for this reduced ordered moment. We have also discussed 
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differences in sum rules between Hubbard and Heisenberg 
models, and refer the reader to more recent work& for ad- 
ditional insights, alternative approaches, and the current 
experimental status. Beyond the single magnon excita- 
tion, we have shown that there is considerable spectral 
weight in the continuum along the zone boundary, and 
that the dynamical structure factor exhibits secondary 
peaks which we understand as arising from the density 
of states of two-spinon excitations. These could be ex- 
plored in experiments with high intensity neutron sources 
in the near future. 
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